#ifndef __BWLIB_NUMBER_THEORY_MILLER_RABIN_H__
#define __BWLIB_NUMBER_THEORY_MILLER_RABIN_H__
#include "../config.h"
#include "module.h"
#include <cstdlib>
#ifdef __BWLIB_WITH_GMPXX
#include <gmpxx.h>
#endif

namespace bwlib
{
	namespace detail
	{
		template<typename NumType>
		inline bool witness(const NumType& p, const NumType& a, 
			const NumType& tq, const NumType& q, NumType k)
		{
			NumType x = mod_pow(a, q, p);
			if(x == 1 || x == tq)
				return false;
			while(k-- != 0)
			{
				x = mod_mul(x, x, p);
				if(x == tq)
					return false;
			}

			return true;
		}
	}

	/* @brief: Miller-Rabin 素性测试
	   @complexity: O((1 + O(1)) log n)
	   @param: p,      待测试数
	           times， 测试次数
			   rnd，   随机数生成器
	   @return: true， p 可能是素数
	   			false，p 不是素数 */
	template<typename NumType, typename RandFactor>
	inline bool miller_rabin(const NumType& num, int times, RandFactor rnd)
	{
		if(num == 1)
			return false;
		if(num == 2)
			return true;
		if((num & 1) == 0)
			return false;
		
		NumType tq = num - 1;
		NumType k = 0, q = tq;
		while((q & 1) == 0)
		{
			q >>= 1;
			++k;
		}

		while(times--)
		{
			NumType a = rnd() % tq + 1;
			if(detail::witness(num, a, tq, q, k))
				return false;
		}

		return true;
	}

	template<typename NumType>
	inline bool miller_rabin(const NumType& num, int times)
	{
		return miller_rabin(num, times, std::rand);
	}

#ifdef __BWLIB_WITH_GMPXX
	template<>
	inline bool miller_rabin(const mpz_class& num, int times)
	{
		return mpz_probab_prime_p(num.get_mpz_t(), times);
	}
#endif
}

#endif
